Importamos las librerías necesarias
import numpy as np
import scipy as sp
from scipy.sparse import diags
import matplotlib.pyplot as plt
import matplotlib.animation as animation
Preparamos la función de animación para poder mostrar el resultado
from matplotlib import rc
rc('animation', html='jshtml')
def plotAnimation(time, space, solution):
fig, ax = plt.subplots()
ax.set_xlabel('x')
ax.set_ylabel('tª')
plotLine, = ax.plot(space, np.zeros(len(space))*np.NaN, 'r-')
plotTitle = ax.set_title("t=0")
ax.set_xlim(np.min(space)-0.1, np.max(space)+0.1)
ax.set_ylim(-1, np.max(solution[0])+0.1)
def animate(t):
pp = solution[t]
plotLine.set_ydata(pp)
plotTitle.set_text(f"t = {time[t]:.1f}")
return [plotLine,plotTitle]
anim = animation.FuncAnimation(fig, func=animate, frames=np.arange(0, len(solution)-1, 1), blit=True)
return anim
Definimos la función usando matrices, aunque luego impondré las condiciones a mano:
def Wave(f, df, Np, Mp, condition, t0=0.0, tn=20.0, x0=0.0, xn=10.0):
N=Np
M=Mp
a=t0
L=tn
#Generamos los nodos
k = (L-a)/N
h = (xn - x0)/M
nodost = [a+k*i for i in np.arange(N)]
nodosx = [x0+h*i for i in np.arange(M)]
#Preparamos la matriz tridiagonal para la derivada segunda en espacio
diagonals = [-2*np.ones(M), np.ones(M-1), np.ones(M-1)]
mat = diags(diagonals, (0, 1, -1)).toarray()
#Aplicamos la condición inicial a los nodos
inicial = np.array([f(x) for x in nodosx])
#Calculamos el primer valor a utilizar en el cálculo usando la derivada
anterior = np.array([f(x) + k*df(x) for x in nodosx])
#Iteramos para calcular las soluciones en cada tiempo a partir de las dos anteriores e imponemos condiciones
result=[anterior, inicial]
for i in np.arange(2, N+1):
actual = 2*result[i-1] - result[i-2] + ((k/h)**2)*mat@result[i-1]
#Condiciones tipo Neumann
if(condition=="Neumann"):
actual[0] = actual[1]
actual[M-1] = actual[M-2]
else:
#Condiciones tipo Dirichlet
actual[0] = actual[M-1] = 0
result.append(actual)
return nodost, nodosx, result
def onda(x):
return np.sin(x*np.pi/10)
def derivada(x):
return np.cos(x*np.pi/10)*np.pi/10
nodost, nodosx, result = Wave(onda, derivada, 500, 100, "Dirichlet")
plotAnimation(nodost, nodosx, result)
nodost, nodosx, result = Wave(onda, derivada, 500, 100, "Neumann")
plotAnimation(nodost, nodosx, result)
He elegido esta condición inicial (a base de ir refinando para ponerla a la izquierda, quizá es algo complicada), a diferencia de la anterior se introduce cierta oscilación extraña, como no parece inestable y no ocurría en el otro caso, puede tratarse de ruido numérico procedente de la aproximación de la condición inicial y su complejidad.
def onda(x):
return 7 * np.exp(-1*x**2)*np.sin(x*np.pi/10)
def derivada(x):
return (-1/10)*np.exp(-1*x**2)*(140*x*np.sin(x*np.pi/10) -np.cos(x*np.pi/10)*np.pi*7)
nodost, nodosx, result = Wave(onda, derivada, 500, 100, "Dirichlet", tn=40)
plotAnimation(nodost, nodosx, result)
nodost, nodosx, result = Wave(onda, derivada, 500, 100, "Neumann")
plotAnimation(nodost, nodosx, result)